
library(mosaic)
library(Stat2Data)
data(SATGPA)

#get the oOriginal correlation

originalr=cor(MathSAT~VerbalSAT,data=SATGPA)
originalr


#Generate lots of correlations where GPAs are randomly permuted
N=10000
randcorr=do(N)*cor(shuffle(MathSAT)~VerbalSAT,data=SATGPA)


#Display of randomization correlations
histogram(randcorr$cor,main="Randomization Correlations")

#See what proportion of randcorrs are beyond the original r
upper=sum(randcorr$cor>abs(originalr))     #count those results > originalr
lower=sum(randcorr$cor<(-abs(originalr)))    #count the lower tail results
pvalue=(upper+lower)/N
pvalue 

